Fluctuations in the relaxation dynamics of mixed chaotic systems 



Roy Ceder and Oded Agam 1 

The Racah Institute of Physics, The Hebrew University, Jerusalem 91904, Israel 

(Dated: October 16, 2012) 

The relaxation dynamics in mixed chaotic systems are believed to decay algebraically with a 
universal decay exponent that emerges from the hierarchical structure of the phase space. Numerical 
studies, however, yield a variety of values for this exponent. In order to reconcile these result we 
consider an ensemble of mixed chaotic systems approximated by rate equations, and analyze the 
fluctuations in the distribution of Poincare recurrence times. Our analysis shows that the behavior 
of these fluctuations, as function of time, implies a very slow convergence of the decay exponent of 
the relaxation. 
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I. INTRODUCTION 

Recurrences play an important role in physics. The 
statistics of the recurrence times of particles govern 
the transport properties of open systems such as quan- 
tum dots attached to external leads, and the relaxation 
characteristics of distribution functions in closed sys- 
tems, see Ref. [l|. Additionally, in the short wavelength 
limit, quantum mechanical properties such as energy 
level correlations^, weak localization- and shot noised, 
are also determined by recurrences of the underlying clas- 
sical dynamics. 

One of the basic statistical characteristics of recur- 
rences in open systems is the distribution of Poincare 
recurrence times, F(t). This is the probability of a tra- 
jectory to return, at time larger than t, to a predefined 
region in phase space. In this paper we focus our at- 
tention on the the distribution of Poincare recurrence 
times of "mixed" chaotic systems (more precisely-one di- 
mensional symplectic maps). In these generic systems 
the phase space consists of islands of regular dynam- 
ics immersed in a sea of chaotic behavior, and around 
these islands there are smaller satellite regular motion 
islands around which are even smaller islands, and so 
forth ad infinitum^. In such systems, a trajectory within 
the chaotic component of the phase space may stick to 
the regular islands for an exceedingly long time. This 
feature is believed to manifest itself in an algebraic de- 
cay of the distribution of Poincare recurrences in the long 
time asymptotic limit: 



F(t) ~ r">". 



(i) 



Many studies have investigated this algebraic decay, 
focusing on two main issues: The question of universal- 
ity of the decay exponent 70, and the calculation of its 
actual value. Yet, after three decades of studies, the an- 
swers to these questions are still controversial. Theoreti- 
cal and numerical studies^— yield a variety of values for 
70 ranging from 1 to 3. Some of the numerical studies 
have obtained different values even for the very same sys- 
tem. A possible explanation to this odd situation is that 
the long time asymptotic behavior of the Poincare recur- 
rence distribution cannot be reached within the existing 



computational power. The numerical calculation of 70 in 
the long time asymptomatic limit is exceedingly difficult. 
Trajectories of a particle moving in the intricate hierar- 
chical structure of the phase space are very sensitive to 
numerical noise. Such a noise may transfer the particle 
from regions of chaotic motion to regular ones and vice 
versa, and can also help in crossing cantori which serve as 
leaky barriers within the phase space. These uncontrolled 
artifacts, apparently, hinder the identification of the true 
asymptotic behavior of F(t). This explanation, however, 
relies on the assumption that even in the absence of nu- 
merical noise the convergence of the Poincare recurrence 
distribution to its asymptotic behavior is very slow. The 
aim of the present work is to study this aspect of the 
relaxation problem. We show that F(t) exhibits time 
dependent fluctuations which do not decay fast enough 
in time. Consequently F(t) at different time intervals 
seems to exhibit different relaxation behavior which may 
be interpreted as different relaxation exponents. This 
behavior is demonstrated in Fig. 1 where the numerical 
calculation of F(t) for the standard map is depicted on a 
log-log scale. It shows that different time intervals may 
be associated with different decay exponents. Thus con- 
vergence to the asymptotic behavior is extremely slow, 
and not monotonous, but rather oscillatory in log t. 

These variations in the relaxation dynamics emerge 
from random-like local variations in the phase space of 
the chaotic component of mixed systems. The trajecto- 
ries that contribute to F(t) as time progresses are tra- 
jectories which become closer to the boundaries of the 
regions with regular motion. The closer the trajectory 
approaches to such a boundary, the longer it sticks to 
it. This behavior implies that self-averaging of the re- 
laxation dynamics is not very effective and therefore the 
convergence to the asymptotic value of the decay expo- 
nent is rather slow as demonstrated in Fig. 1. 

Obviously, a direct study of the fluctuations of F(t) 
around its asymptotic behavior ([T]), whether analytical 
or numerical, suffers from the same difficulties of the cal- 
culation of the asymptotic decay exponent itself. To cir- 
cumvent this difficulty it is convenient to introduce an 
ensemble of mixed chaotic systems and to use ensem- 
ble averaging in order to extract the properties which 
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FIG. 1: A log-log plot of the survival probability of the stan- 
dard map showing local fluctuations in its decay power 7. 

characterize the relaxation dynamics. This is similar in 
spirit to disorder averaging— which allows one to char- 
acterize the dynamics of a particle in disordered system, 
such as the correlation function of its position at dif- 
ferent times. However, unlike disordered systems where 
the meaning of their ensemble is well understood, it is 
unclear what might be the invariant measure of the en- 
semble of mixed chaotic systems. Nonetheless, the role 
of the ensemble averaging which we shall employ here 
is merely a regularization procedure which allows one to 
extract the intrinsic properties of the pure system. It is 
similar to the introduction of an infinitesimal noise into 
the dynamics of a hard chaotic system in order to ex- 
tract the Pollicott-Ruelle resonances which characterize 
the relaxation of the pure system^. Thus here we shall 
use disorder diagrammatics in order to perform averages 
and identify the correlation function of the fluctuations 
in the return probability. Then we shall use this result 
to calculate the typical value of the local fluctuations in 
the decay exponent S"f = 7 — 70 , and show that 

(<57 2 )~i- 2 7(logi) (2) 

where e <C 1, and /(logt) is an oscillatory function. 

Our analytical analysis of the problem is based on the 
rate equation model for the dynamics of mixed chaotic 
systems developed by Meiss and Otl£. In Sec. II 
we present the model where following Cristadoro and 
Ketzmericki^ we add a small random component to the 
transition rates of the original Meiss-Ott model. In sec- 
tion III we present the solution of the pure model, and in 
section IV we use the results of disorder diagrammatics 
to calculate the correlation function of the fluctuations in 
the return probability. This correlation function will be 
used in Sec. V in order to derive formula @ for the local 
fluctuations in the decay exponent. In Sec. VI we shall 
present numerical calculations which support our ana- 
lytical results, and conclude in Sec. VII. The technical 
details of our calculations can be found in the Appen- 
dices. 




FIG. 2: The hierarchical structure of the phase space of mixed 
chaotic systems (left panel), and the Meiss-Ott tree model 
(right panel). 



II. THE RANDOM TREE MODEL 

In this section we present our model for the ensem- 
ble of mixed chaotic systems. It is based on the as- 
sumption that the statistical characteristics of the dy- 
namics of mixed chaotic systems can be described by 
rate equations. First, let us recall the rate equations 
approach for mixed chaotic systems introduced by Meiss 
and Ott-. Their model assumes that the dynamics within 
the chaotic component of the phase space can be approx- 
imated by the Master equation for the probabilities of 
finding the particle in states which respect the self-similar 
structure of the phase space, see illustration on the left 
panel of Fig. 2. Each one of the Meiss-Ott states is as- 
sociated with a definite region of the phase space whose 
boundaries (represented by the dashed lines in the left 
panel of Fig. 2) are determined by lowest flux cantori 
encircling the relevant set of regular islands. 

The topology of the phase space division is that of a 
Cayley tree, as shown in the right hand side of Fig. 2. 
Namely, a particle in a given state might move either 
to a single lower state of the hierarchy or to one of two 
possible higher states. A "level" transition is associated 
with transition to a state which is closer to the boundary 
circle of the island chain it is presently revolving around, 
while a "class" transition corresponds to the case where 
the particle moves into a higher state associated with one 
of its present satellite island chains, see Fig. 2. 

The binary structure of Meiss-Ott tree allows one to 
designate the states of the system by binary numbers: A 
state reached by level transition is denoted by the number 
of the previous state to which a figure " 1" is added, while 
a that reached by class transition is obtained by adding 
the figure " 0" , as demonstrated in the right hand side of 
Fig. 2. Let n represent some arbitrary state on the tree, 
and denote by Dn the nearest state at the upper part 
of the tree from which the particle may arrive. We shall 
also denote by nO and nl the two states down the tree 
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reached by class and level transitions, respectively. With 
these definitions, the Master equation takes the form 



dP, 
~dt 



(W, 



W n ^ nl )P n 



+ W n0 ^ n PnO + W nl ^ n P nl + W D P D n, (3) 



where P n is the probability to find the particle in the state 
n, while W n ^ m is the transition rate from the state n to 
to. Meiss and Ott assumed that these transition rates 
satisfy a simple scaling behavior: 



— £ 777 — £l 



W Dr . 



and 
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W n ^nO So W n ^ n \ 



El 



(4) 



(5) 



where eo, ex, ujq, and u>±, are constants which have been 
estimated to be 8 -: 



e ~ 0.143, e x ~ 0.382, 
w ~ 0.0142, ujx ~ 0.0532, 



(6) 



Within this model it is also assumed that a particle stay- 
ing in the upper state of the tree leaves the tree at rate 
r and never return back. From here on, the escape rate 
r will be set to unity by choosing proper units of time. 

It is known, however, that the rates ratios |@J and ((SJ) 
fluctuate considerably at different positions of the tree^i. 
Therefore it is natural to generalize the Meiss-Ott model 
by adding a random component to the transition rates, 
i.e. to replace the constant transition rates by fluctuating 
ones: 



W n 



I (1 + Cnm) 



(7) 



where £ nm are assumed to be uncorrelated random vari- 
ables with zero mean, 



(8) 



Here a is a dimensionless constant which controls the 
amount of randomness of the ensemble. 

This ensemble is based on two main assumptions: (a) 
The main source of fluctuations comes from the flux ex- 
change between states, while deviations from the scal- 
ing behavior of phase space areas of the states may be 
neglected, (b) The fluctuations in the flux exchange 
through different boundaries of the states are uncorre- 
lated. The above choice implies that statistical proper- 
ties of the dynamics on small scales of the phase space 
is the same as at larger scales with the proper rescaling 
of time. We shall refer to this ensemble as the "random 
tree model" . 

In order to identify the relation between the functions 
P n (t) and the survival probability, F(t), let us sum the 
Master equation ^ over all the states of the tree n: 



dt 



dF(t) 
dt 



Px 



(9) 



The left hand side of this equation is precisely the deriva- 
tive of the survival probability, F(t) = ^ n P n , since the 
sum over P n is the probability to find the particle at any 
site on the tree. From here and (Q]) it follows that the 
long time asymptotic decay of Px (t) is 



Pi(t)~f 



-I-70 



(10) 



Thus one may study the convergence to this asymptotic 
limit by analyzing the behavior of the return probability 



P(t) = Px(t) where P„(0) = 5 n ,x 



(11) 



This is the probability density that a particle is found at 
upper site of the tree, n = 1 assuming that at t = the 
particle is placed at the same site. 



III. SOLUTION OF THE NONRANDOM 
MODEL 

In this section we review the solution of the Meiss-Ott 
model that describes the pure system, thereby presenting 
some of the ingredients which will serve us in the next 
section where we consider the random model. Let us de- 
fine the probability P n ,m{t) of a particle, initially placed 
at site m, to be at site n at time t. We define the Green 
function, G„, m (s), as the Laplace transform of P n ,m(t)i 



dte- st P n , m (t). 



(12) 



Our aim is to characterize the analytic structure of the 
green function associated with the return probability to 
the upper site of the tree, Gx,i{s). In particular we 
shall focus our attention on the analytic structure in the 
vicinity of the point s = which governs the long time 
asymptotic behavior of the return probability (jlip . In 
this vicinity one expects Gx,x(s) to take the form 



Gx,x(s) = a(s)+b(s)s^ + 



(13) 



where a(s) and b(s) are some analytic functions. Since 
the inverse Laplace transform of an analytic function de- 
cays, in time, faster than any power law, the long time 
asymptotic behavior (|10p comes from the nonanalytic 
contribution represented by the second term of the ex- 
pansion flT3|) . Meiss and Ott proved that 7 satisfies the 
dispersion equation 8 : 



U! £ a 7 + lu 1 e 1 1 



1. 



(14) 



This equation is obtained by substituting (1131) into the 
equation that the Green function Gx,x(s) satisfies. It is 
rederived using a diagrammatic approach in Appendix 
A. 

Let us examine the solutions of the dispersion equation 
(TT41 . Apart from a single purely real solution the other 
solutions of this equation appear in complex conjugate 
pairs. In Fig. 3 we present some of these solutions in the 
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FIG. 3: The solutions of the dispersion equation (fTl)) 

complex plain (Re(7), Im(7)). From here it follows that 
the solutions are divided into two groups. In each one of 
these groups, the real parts of the solutions are approxi- 
mately the same. The difference between the solutions is 
mainely due to their imaginary parts. The solution with 
smallest real part is 

7o = 1.96424 (15) 

This solution governs the asymptotic long time behavior. 
The solutions with the next smallest real part are 71 = 
1.96469 ± z6. 47555, while the next solutions are 72 = 
1.96605 ± il2.9511i. Thus the general expansion of the 
Green function near s = is 

Gi,i(«) = a{s) +J2K{s)s^ + ■ ■ ■ (16) 
v 

The form of decay associated with complex solutions of 
the form s lR±ZJI is t~ lR cos(7/ log t + <f>) where (j) is some 
constant. Thus neglecting the difference in the real parts 
of the solutions we see that the long asymptotic behavior 
of the return probability of the pure system is 

P{t) = P l {t)~t->°- 1 h l {\ogt) (17) 

where h\(x) is non-periodic oscillatory function whose 
behavior depends on the initial conditions of the system. 
hi(x) may be approximated by the leading solutions of 
the dispersion equation (fl"4"|) , i.e. 

hi(x) w ao + a± cos (7 a: + <fi) (18) 



where ao, a\ and <f> are positive constants, and 
7' = 6.47555 is the imaginary part of the next to leading 
solution of the dispersion equation (fl"4f . 



IV. SOLUTION OF THE RANDOM TREE 
MODEL 

While ensemble averaging is not required for extracting 
the long time asymptotic decay of the return probability, 
the sample to sample fluctuations are defined only with 
regard to ensemble. In this section our goal is to cal- 
culate the correlation function of the sample to sample 
fluctuations: 

C(t,t') = (6P(t)5P(t')) (19) 

where 

SP(t) = P(t) (P(t)) (20) 

is the deviation of the return probability of the system 
with the random component, P(t), from the average over 
the ensemble (P(t)). 

Since we employ the ensemble averaging procedure 
only as a tool for extracting the intrinsic properties of 
the system, we may assume that fluctuations in the tran- 
sition rates are small, a <C 1, and exploit a as the small 
parameter of the perturbation theory. The details of this 
perturbation theory are presented in Appendix B. 

The correlation function (|19l) may be written as a dou- 
ble inverse Laplace transform 

C(t,t')= / / ^Le^QM (21) 

where the constant c in the integration contour is set to 
ensure the convergence of the integral, and Q(s', s) is the 
disconnected part of the correlation of the Green function 

Q(s',s) - (G 1A (s')G 1A (s)) - (G 1A (s'))(G hl (s)), (22) 

where Gi i i(s') denotes the Green function of the random 
tree model. In appendix B it is shown that, to the lowest 
order in a, Q(s's) stratifies the equation: 



Q(s',s) = A(s',s) + Gl 1 (s')Gl 1 (s) 



(23) 



where 



and ctj (s', s) is a function expressed in terms of the Green 



A( S ', s ) = a 2 G? a ( s ')G? a ( S ) 



1 + w l a i( s '' s ) 
j=o,i 



(24) 
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function of the pure system, G\.\ (s). Its explicit form can 
be found in Appendix B. 

From the structure of Eq. (|23|) it follows that its so- 
lution contains two contributions: The homogeneous so- 
lution and the inhomogeneous one. Consider first the 



inhomogeneous solution which we denoted by Qi n {s's). 
Focusing our attention on the terms which are relevant 
for the long time asymptotic behavior we expand this 
solution in the form 



Qin(s', a) = g<°> (s\ s) + J2 (<lP W, s)^- 1 + e ] is, s')^- 1 ) + £ g« (a', S )s'^^s^ + ■■■ (25) 



where (s', s) represent functions that are analytic in s 
and s', and 7„ are the solutions of Eq. (TT4")) . Substituting 
this form in (|23|) it is evident that with a proper choice of 
the functions (s' , s) one can satisfy the equation term 
by term. The first two terms of the expansion (1251) are an- 
alytic either in s or in s' and therefore do not contribute 
to the long time asymptotic behavior of the fluctuations. 
The third term, near a = s' = 0, may be written in the 

form of a product: x^ 2 (E„ M ) 57 "" 1 )^ M0)s' 7 " -1 ) 
where b„(0) are the expansion coefficients of the Green 
function (|16[) . and x is a constant which may be expressed 
in terms of a(0), and the parameters ©. This implies 
that the contribution to the correlation function of the 
fluctuations in the long time asymptotic limit is propor- 
tional to the product of the return probabilities: 

C in (t,t') = X v 2 P(t)P(t') (26) 



where P(t) is the return probability of the system with- 
out the random component, see Appendix A. 

Consider now the homogeneous solutions of Eq. (|2"3"|) , 
which we shall denote by Qh(s's). The total solution 
of the equation is Q(s',s) = Qi n (s's) + Qh(s's) and its 
inverse Laplace transform gives the correlation function 
C(t, t')-see Eq. (l2~Tj) . The homogeneous solution is needed 
in order to satisfy the condition 

C(0,t') = £0,0) =0 (27) 

for any t and t' . This is because at the initial time of 
evolution 5P(0) = as the system is prepared in such a 
way that the particle is with probability one at the upper 
state of the tree, "1". To obtain the long time asymp- 
totic behavior of the homogeneous solution, we expand 
Qh{s' , s) in the form 



Q h (s\ a) = w<°> (a', a) + (vP (s',s)st + (a, a')s>») + £ (*', s)a^ a'^ +■■■ (28) 



where v^\s' , s) are analytic functions at a = a' = and 
fi u are unknown exponents. Substituting this expansion 
in the homogeneous equation obtained from (|23l) . and 
solving the resulting equation term by term, one obtains 
that the contribution associated with the slowest decay 
exponents (which comes from the third term in the ex- 
pansion (|28|) should satisfy the dispersion equation: 

(woe^) 2 + {^e^) 2 = 1 (29) 

where 



This equation, similar to Eq. (Q3J) , has many solutions for 
jl which form two branches similar to those which appear 
in Fig. 3. The solution with the smallest positive real part 
satisfies the relation /2 > 70, and for the parameters © 
its value is jlo = 2.13852. Other exponents belonging 



to the same branch have the approximate form p,o ± ifj/j 
for example n[ = 3.23377 and fi' 2 = 6.46753. Since the 
dispersion equation (1291) constrains only the sum of ex- 
ponents (|30p there is apparently a continuous set of so- 
lutions associated with different values of the difference 

between the exponents, /i„ — [i v < , the weight of each one 

(2) 

of these solutions is determined by v v v , (s',s). This ap- 
proach does not allow us to determine these weights since 
our solution is valid in the long time asymptotic limit 
while the condition (|2"T|) corresponds to short time dy- 
namics. Nevertheless, we may express the homogeneous 
contribution to the correlation function, in the long time 
asymptotic limit, as a sum 

C h (t,t') = Rc^cw, ( t -V 2 ^-" + ft-vpi^-v) (31) 

where c v ^ are some unknown constants and \x v are the 
solution of the dispersion equation (|2"9"|) . Here we assume 
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that the sum is only over the solutions with positive imag- 
inary part. In particular for t' = t one obtains 



where dj and <pj are some constants, and represent the 
imaginary part of the solutions of the dispersion equation 
(EH). 

Now from (|26|) . (|32|) and the results presented in Sec. 
V we get that the variance of the normalized fluctuations 
Sp(t) = SP(t)/(P(t)), in the long time asymptotic limit, 
and to the leading order in a (where (P(t)) should be 
replaced by P)takes the form 

(Sp 2 (t))= x <r 2 (l + f3i(t)) (34) 



where 



r m ^ h2 
x K 



e = p - 7o ^ 0.174 (36) 



Thus the normalized fluctuations in the return proba- 
bility decay to a constant value very slowly and in an 
oscillatory manner as function of log t. Finally, using the 
above results for the homogeneous and inhomogeneous 
contribution to the correlator C(t + At/2), t — At/2), and 
expanding it to second order in At/t, we obtain: 



C h (t,t) = (SP(t) 2 ) h ~ a 2 r 2 ^h 2 (\ogt) (32) 

where h 2 (log t) is an oscillatory function which similarly 
to h%(t) (see Eq. (ITS)) ~) may be approximated in the long 
time asymptotic limit as 

and 

h 2 (x) « d + di cos(2fi' 1 x + 0i) + d 3 cos(2fi' 2 x + </>) + • -(33) 



„ ( At At 



xt 



-270 



h 2 (\0gt) 



(37) 



r 



where 

hi(x) = la h 2 1 (x)-h' 2 (x) + h 1 (x)(h';(x)-h'(x)) (38) 

and is h 2 (log t) is an oscillatory function with expansion 
similar to (|33|) but with different coefficients. From here 
one obtains the correlator of the normalized fluctuations 



(Sp(t 



2 I 



Spit- 



2 



)) 



where 



(SP 2 (*)> 



/3 2 (t)=t- 



( At 



i - - 2 m { t 



; 7o/i 2 (logf) - h 2 (logt) 
2 X h 2 (\ogt) 



(39) 



(40) 



To obtain this result we take into account that hi(x) is 
approximately constant and thus h\{x) may be approxi- 
mated by the first term on the right hand side of 



where (P(t)) is the long time asymptotic limit of the av- 
erage of the return probability over the ensemble which 
to the leading order in the disorder may be taken as the 
return probability of the pure system, (|17l) . The function 
8p(t) = P(t) / (P(t)) — 1 is the normalized sample to sam- 
ple fluctuation which is characterized by the correlation 
function (|3"9")l . Assume one wishes to extract the decay 
exponent 7, of some particular system, from the mea- 
surement of P(t) at two time points, say t ± At/2 where 
the time difference At between the measuring points is 
assumed to be smaller than t, At <C t (as is the general 
situation when one tries to extract the limiting value of 
the decay exponent from the tail of P(t)) . Assuming that 
within this range the local behavior is P(t) cx t _1 ~ 7 , one 
obtains that 



t /l + 5p(t- At/2)\ , s 

^ 1 + ^ + At h i l + 6p\t + At/2) )> ^ 



V. SAMPLE TO SAMPLE FLUCTUATIONS OF 
THE DECAY EXPONENT 

In this section we calculate the typical value of the 
sample to sample fluctuation in the local behavior the 
decay exponent 7. For this purpose let us present the 
return probability in the form 



P(t) = (P(t))(l + Sp(t)), 



(41) 



Thus Sj — 7 — 70 can be associated with the normalized 
fluctuations of the return probability, 6p(t). Expanding 
the above logarithm to the leading order in 5p(t), one 
may express the mean square of the fluctuations in the 
decay exponent in the form 



(43) 
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Using (|39j) and (jM) we obtain 

(^}~ X( 7 2 (ft(i)+fX(t)) (44) 

where /?i (£) and /32 (i) are given by (l3"5"j) and (14*01 respec- 
tively. Both functions are oscillatory in log t and decay as 
t~ 2e . Thus the typical fluctuation in the value of decay 
exponent decays as £7 ~ t~ e , where e « 0.174. From here 
one concludes that the fluctuations decay very slowly in 
time and in an oscillatory manner. 

VI. NUMERICAL STUDY 

In this section we present our numerical study which 
has the following goals: First to verify our analytical so- 
lution. Second to show that the results that we have 
obtain in the limit a <C 1 are insensitive to the magni- 
tude of the random component of the transition rates. 
Third, to compare our results with those obtained from 
the exact dynamics of an ensemble of symplectic maps. 

Our numerical study of the random tree model is per- 
formed directly for the survival probability, F(t), rather 
than the return probability P(t). The computations are 
performed using 1000 realizations of the ensemble of ran- 
dom tree model with 9 generations. The random vari- 
ables £ n ,m, are chosen form a uniform distribution. 

The relation (j9)) implies that F(t) behaves in a similar 
way to P(t). In particular one expects that 

F(t) = r 7D /i(logt) (45) 

where similar to (|T7|) . the function h (log t) is an oscil- 
latory function whose expansion is the same as that of 
hi(\ogt) but with different coefficients. Yet, formula (fT"5)) 
describes the behavior in the very long time asymptotic 
behavior, while within the time limits were our calcula- 
tion is reliable one has to take into account also terms 
associated with the second branch of solutions of the dis- 
persion equations (TT41 that are associated with faster de- 
caying rates (see right branch of solutions in Fig. 3). 

The numerical calculation of h(\ogt) for the pure tree 
model is presented by the black dots in Fig. 4. This solu- 
tion is obtained by calculating the normal modes of the 
tree model and expanding the solution in these modes. 
The solid line in this figure represents the analytical so- 
lution described by the function 

h(x)~a + — + a 2 t 10 ' 1 ' 1 cost's + fa) 

+ a 3 f ^ cost's + fa) (46) 

where j[ ± vy'{ = 1.9646 ± i6.4755 and i 2 ± = 
2.4097 ± i3.223 are the slowest oscillatory solutions of 
the dispersion equation (fT4"f while the term a\jt comes 
from the linear in s term of the expansion of the function 
a(s) of the Green function (|13p . The parameters a,j and 
fa are fitting parameters. The analytical solution gives 
an excellent fit to the numerical data except at very large 



times where numerical errors are large. Notice that the 
amplitude of the oscillatory component of h{t) is very 
small, thus to good approximation it may be considered 
to be constant. 



« 102 



1.01 




10 2 10 3 10 4 10 : 



t 

FIG. 4: The numerical (black dots) and analytical (solid line) 
results for the function h(t) which dictates the oscillatory na- 
ture of the survival probability (|45[) . 

Let us consider the behavior of (Sj 2 ). From formula 
(|4*4"| . expressions (1331) and (|4"0"1) . and taking into account 
that to a good approximation h(t) may be replaced by 
constant, it follows that t 2e (Sj 2 ) may be expanded in 
terms of the solutions \J ± , of the dispersion equation 
(JUJ), namely 

t 2t (Sj 2 ) = t 2fM> ^2 v 3 t ' t "' 3 cos(/i"j log* - fa) (47) 

3 

The black dots in Fig. 5 represent the numerical calcu- 
lation of this function for the case were a = 0.04, while 
the solid line is the analytical expression (|4T|) in which 
we took the four slowest oscillatory solutions of (|2"9")l and 
use Vj and fa as fitting parameters. The very good agrc- 
ment between the numerical and analytical results shown 
in Fig. 5 proves that the typical amplitude of the fluctu- 
ation, #7 decays as t~ e where e is given by (f5B|) . It also 
shows that the oscillatory behavior of (Sj 2 ) is dictated 
by the oscillatory solutions of (|29|) . 

Our analytical results have been derived for weak disor- 
der. In this limit the strength of the random component 
in the transition rates appears only through the prefac- 
tor which controls the magnitude of the fluctuation (see 
Eq. (|4"4"|) ), while the functional dependence on time is 
independent of the a and depends only on the intrinsic 
properties of the pure Meiss-Ott model. Our numerical 
comparison, shown in Fig. 5, has been also performed for 
weak disorder <j = 0.04. It is, however, worthwhile to 
clarify to what extent these results depend on the value 
of cr. . To this end we compute (<57 2 ) for various values 
of cr, between 0.04 and 0.2 (Notice that the higher value 
is not far from the upper limit of the widest possible 
uniform distribution for which a = l/v!2 ~ 0.29). 
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FIG. 5: The numerical (black dots) and analyticl (solid line) 
results for the oscillatory factor of the sample to sample fluc- 
tuations of the decay exponent (|4T[) . 

The functions t 2e (8^f 2 ) / a 2 are presented in Fig. 6. The 
fact that they almost collapse on the same graph, implies 
that our results are almost independent of the strength 
of the random component of the transition rates. 
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FIG. 6: (Color online) The functional behavior of the sample 
to sample fluctuations of the random tree model for various 
values of the disorder strength. 

Finally it is instructive to compare our results with 
those of the true dynamics of mixed chaotic systems. 
For this purpose we study an ensemble of chaotic maps 
obtained by adding a small random component to the 
standard map: 

In+i = I n + Ksm(9 n )+K{9 n ) 

n +i = n + In+i mod(27r) (48) 

Here K is the kicking strength, and 11(9) is a Gaussian 
random function with zero mean and periodic correlation 

(11(9)11(9')) = a 2 ^ e- ie -°'T m) * (49) 

m 

The constant a controls the strength of the random con- 
tribution, while I is the correlation length of 11(9). Thus 
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FIG. 7: (Color online) The relative fluctuations in the distri- 
bution of Poincare recurrences of the random standard map 
model ((48} versus time. 

each realization of 11(9) corresponds to a slightly differ- 
ent symplectic map, which may be viewed as a different 
member of the random tree model. Thus averaging over 
this ensemble is expected to have similar effects to aver- 
aging over the transition rates of the random tree model. 
For the numerical computation we choose a = 0.21 • 1CP 5 
which is much smaller than the largest value of K, for 
which the unstable fixed point at / = and 9 — n re- 
mains unstable. The initial points of the various trajec- 
tories of the particles, are chosen to be in the vicinity 
of this fixed point (within a square of size 10~ 6 x 10 -6 ), 
and the particle is assumed to leave the system when 
it crosses the line I — 0. We choose kicking strength 
K = 0.971635406, correlation length I = 0.2, and the cal- 
culation is performed using a quadruple precision. The 
map was iterated up to 5 • 10 5 iterations, which is the 
range within which the results agree with those obtained 
by double precision. For each realization of the random 
function 11(9), the survival probability was calculated us- 
ing 2- 10 9 trajectories. The correlations of the fluctuation 
have been obtained by averaging over 260 different real- 
izations of the random component. 

In Fig. 7 we present the results for the fluctuations 
in the decay exponent, (S^ 2 ), as function of \ogt. This 
graph shows an oscillatory behavior of the fluctuations 
but with no apparent decay, namely e 0. 

VII. CONCLUSIONS 

In this work we have extended the rate equation ap- 
proach of Meiss-and Ott in order to calculate the func- 
tional behavior of the sample to sample fluctuations in 
the local decay exponent 7 of the survival probability. 
Our calculations were conducted within the leading or- 
der perturbation theory in the strength of the random 
component. In this limit, it has been shown that apart 
from the amplitude of the fluctuations in 7, their func- 
tional dependence on time is dictated by the intrinsic 
parameters of the Meiss-Ott model. Moreover, our nu- 
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merical study shows that this conclusion is valid also at 
strong disorder. Our approach cannot give the typical 
value of the random component, a. However, numerical 
evidence from the study of mixed chaotic systems sug- 
gests that it is of order unity. If we also assume that the 
typical magnitude of the fluctuations in 7 is character- 
ized by the decay exponent e = 0.174, then even after a 
million iterations the typical fluctuation in 7 is of order 
of one tenth. 

Moreover, our numerical study of the ensemble of sym- 
plectic maps (j48|) . does not show any sign of decay of the 
fluctuations. Namely, within the limits of reliable nu- 
merical calculation it is difficult to obtain the asymptotic 
value of the decay exponent, 70 (as also evident from Fig. 
1). One may speculate that this is because the assump- 
tion that the random component in the transition rates is 
uncorrelated is not correct, and that its effect cannot be 
considered to be perturbative. For instance it is plausible 
that strong and correlated random components may con- 
fine the dynamics, in the long time limit, to a very small 
number of branches of the tree model. This implies that 
self averaging is not effective and therefore fluctuations 
in the normalized return probability, 8p(t) do not decay 
in time, or in other words e ~ 0. This result is also ob- 
tained for a very asymmetric random tree model where 
class transitions are much smaller than level transition 
Eq/sijLUo/uji — > 0. In this very asymmetric model, e — > 
and therefore does not decay in time. 

Let us, finally, remark that the tree model of Meiss 
and Ott represents an uncontrolled phenomenological ap- 
proximation of the exact dynamics which is based on the 
assumption that within the phase space region associated 
with a given state, the relaxation is much faster compared 
to the transition time to other states. The status of va- 
lidity of this assumption, however, is unclear—. Never- 
theless, our analysis corresponds to a worst case scenario 
were self averaging is effective, and even in this case we 
obtain that the fluctuations in return probability decay 
extremely slowly in time. Thus the convergence of the 
survival probability, F(t), to its asymptotic value is very 
slow. This explains the wide range of results for the de- 
cay exponent 7 obtained by numerical studies [5-18] in 
the last three decades. 
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VIII. APPENDIX A: DERIVATION OF EQ. (Tl4)l 

In this appendix we derive the dispersion equation (|14j) 
which determines 70 as well as the other exponents of 
the relaxation modes shown in Fig. 3. For this purpose 
we first derive an exact equation of the Green function 
Gi^i(s) and then show that a solution of the form (fT5|) 
leads to the dispersion equation. 




FIG. 8: Diagramatic expansion of the Green function of the 
Meiss-ott tree model 

From the Master equation ([3]) and from the definition 
of the Green function (fT2"j) we obtain: 

((3 + Hn)Sj. n ~ \W n ->-j) G n . m (s) = Sj <m (50) 

n 

where 

k 

and A is a dummy parameter introduced for the expan- 
sion of the Green function as a perturbation series, and 
should be set to one at the end of the calculation. In 
particular to the zeroth order in A, G n ,m{s) ~ gn(s)5 n , m 
where 



Now let us expand the Green function G, lj)n (s) to all 
orders in A. The diagrammatic representation of this 
expansion is presented in Fig. 8. Here the thick line rep- 
resents the exact Green function, thin lines represent the 
zeroth order Green function, g n ,m{s), and wiggly lines 
represent the product W n ^ m W m ^ n - Notice that in this 
expansion wiggly lines cannot cross because there are no 
loops on Cayley trees. Thus one may write an equation 
for the Green function Gi,i(s) : 



GiA s ) = 9i(s) + gi{s) ^-►y^SJwWWy-^i.iOO ( 53 ) 

3=0,1 
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where the sum over j is the over the nearest neighbors sites 10 and 11, and Gy^s) denotes the Green function on 
site lj assuming that if the particle reaches site 1 it disappears. From the self similarity of the Cayley tree it follows 
that y(s) and Gx,x(s) are related by a simple rescaling of time: 



G (1) 



W x - 



-Gi, 



W] 



— — Gi i 



Thus the Green function satisfies the equation: 

+ 1 + ujq + vx - ^oGi,i 



wiGi 



(54) 



(55) 



Now substituting (|13[) into this equation gives in a 
straightforward manner equation for a(0), 

[1 + cj + wi - w o a(0) - wio(0)] a(0) = 1 (56) 

whose solution is 

a(0) = 1. (57) 

Now substituting Gx,x{s) = 1 + 6(0)s 7 + ■ ■ • one obtains 
the equation 



6(0) 



s 7 - UJQ 



- wi 







(58) 



which is equivalent to the dispersion equation (|14[) . 



IX. APPENDIX B: DISORDER 
DIAGRAMMATICS 

In this appendix we compute the correlation function 
of two Green functions Gj t k(s). Our main focus is the 



J 



r 



disconnected part of the correlator [221 an d we limit our 
considerations to the the lowest order perturbation the- 
ory in a, for which the result is independent of the precise 
distribution of the random variables 

Taking into account the random component of the 
transition rates we may write the equation for the Green 
function in the form: 



(pnsn( s ) + "«,m) G m> fc(s) 



(59) 



where 



in •// s/n 



(60) 



is the addition to Eq. (|50p which comes from the random 
component of the transition rates. 

From the definition of £ nm (see Eq. ©) we have 
(Sn,m) = and 



n,k^k,m^n, 



r 



(61) 



Notice also that this expression satisfies the relation: 

(Hn,feHi, m ) = (62) 



which follows immediately from the definition (|60j) . 

Consider now the disconnected part of two Green func- 
tions (|2"2"|) . The general diagrams describing the leading 
order expansion of the average (Gx,x(s')Gx,x(s)) (up to 
a 2 ) are shown in the upper panel of Fig. 9. Here solid 
lines represent the exact Green function of the system 
without the random component, while dashed lines stand 



for pairs of random components of the transition rates 
given by ([6"T]) . 

The first tree diagrams shown in panel (a) represent the 
disconnected part of the correlator (Gi ) i)((s')Gi ! i(s)), 
which is of no interest here and therefore to the leading 
order in cr, the connected part of the correlator 



Q(s's) = (SGx,x{s')SGx,x(s)) 



(63) 



is given by the diagram shown in panel (b). The expan- 
sion on the right hand side of this panel refers to the point 
along the tree in which the random contribution is taken 
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(a) 



(b) 



(c) 



1 1 1 n m 1 1 k I 1 

1 n m 1 11 v n in v 1 1 11 v V q n m R v V 1 1 

+ | i' 1 1 ~ + jjj 1 1 + ■ 

i t^i i i miii i nvn m>>n 

1 1 1 1 1/ u 1 i 

~T~ + E [ j j 

1 l ^~ 0,1 111/ LI 1 1 

A(s,'i) 



into account. The first diagram which will be denoted 
henceforth by A(s's) corresponds to the case where the 
random contribution is taken into account along a bond 
from site " 1" to one of its nearest neighbors sites. The 
next diagram represents the case where disorder is taken 
between site j and its nearest neighbor sites, and so on. 
Panel (c) of this figure shows the resummation of all these 
diagrams. From here it follows that 



FIG. 9: Diagrammatic expansion of the correlator 
(Gi,i(s')Gi,i(s)). The dashed line crossed by a segment 
stands for a sum over all diagrams containing one or more 
inner indices equal to the index which appears near the seg- 
ment. 



i iiii y- o, i 
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11 1/ l 
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FIG. 10: The diagram A(s's) 



(SG 1A (s')6G ltl (s))=A(s',s)+ 53 Si4(*) 2 ^i4(^) a Wi-.yW^-.i(* G W,w( a ')*Gy,y(a)) (64) 

3=0,1 

Within the leading order in a one may assume that the Green functions satisfy the same scaling relation (|54]l and 
therefore 




Substituting this relation to (l64l and using definition (|63|) o ne obtains Eq. (l23l). 

Finally let us calculate the diagram A(s',s) whose ex- while dashed lines represent the disorder. Taking into 
plicit expansion is shown in Fig. 10. Here solid lines account the relations: 
stand for the exact Green function of the pure system, 

I 

Gi,ij(s) - G 1 , 1 (s)W 1 ^ lj G^ lj (s) = ^G M ( a )G M (-) (66) 

£ j \ £ j / 

Gii.i(s) = G^^W^G^is) = <2i,i(a)G M (^j (67) 

I 

and using formula (|6ip in order to evaluate the contribu- tion which comes from the the disorder line, leads to Eq. 
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(1221) with 



Gi.i - +G1.1 - 



26',. i ( ^) Gi,i f - 



(68) 
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